# library(devtools)
# install_github('AnneChao/iNEXT.3D')
## import packages
library(iNEXT.3D)
library(readxl)
library(reshape2)
library(doBy)
library(dplyr)
library(plyr)
library(stringr)
library(ggplot2)
library(vegan)
library(ggdendro)
library(patchwork)

ggplot theme setting

large <- theme(legend.title = element_text(size=15),
        legend.text = element_text(size=15),
        axis.title = element_text(size=15),
        axis.text = element_text(size=15),
        strip.text = element_text(size=15))

rotate <- theme(axis.text.x = element_text(angle = 90, hjust=1))

no_strip <- theme(strip.background = element_rect(colour=NA, fill=NA),
                  strip.text = element_text(colour=NA))

Preparing tanaid data

ta0 <- read_excel("../Data/Tanaidacea community.xlsx", sheet = 3)
ta <- ta0[, 1:2]
# Create factor "Microhabitat"
ta$Microhabitat <- str_replace(ta0$Microhabitat, " \\s*\\([^\\)]+\\)", "")
# Create factor "Replicate"
ta$Replicate <- str_extract(ta0$Microhabitat, "(?<=\\()\\d+(?=\\))")
ta <- cbind(ta, ta0[,-1:-3])
# Remove zero rows
ta <- ta[rowSums(ta0[,-1:-3])!=0, ]
ta$Site <- factor(ta$Site, levels=c("Shitiping", "Jihuei", "Jialulan"), labels = c("A", "B", "C"))
ta$Season <- factor(ta$Season, levels = c("SP", "SU"), labels = c("4", "8"))
names(ta)[2] <- "Month"
# Only keep data for analysis
b <- ta[, -1:-4]
rownames(b) <- with(ta, paste(Site, Month, Microhabitat, Replicate, sep="_"))

Prepare CWB data

library(rjson)
library(dplyr)

cwb <-  fromJSON(file="../Data/C-B0050-001.json")

# Meta data 
meta <- lapply(cwb$cwbdata$Resources$Resource$Data$SeaSurfaceObs$Location, FUN=function(x)x[[1]])

# Monthly temperature
dat <- lapply(cwb$cwbdata$Resources$Resource$Data$SeaSurfaceObs$Location, FUN=function(x)x[[2]]$Monthly)

MonthlyTemperatureFun <- function(i){
  cbind(StationNameEN = lapply(meta, FUN=function(x)x$StationNameEN)[i] %>% as.character %>% as.factor,
        cbind(
          StationLatitude=lapply(meta, FUN=function(x)x$StationLatitude)[i] %>% as.numeric,
          StationLongitude=lapply(meta, FUN=function(x)x$StationLongitude)[i] %>% as.numeric,
          DataMonth=lapply(dat[[i]], FUN=function(x)x$DataMonth) %>% unlist %>% as.numeric,
          Maximum=lapply(dat[[i]], FUN=function(x)x$SeaTemperature$Maximum) %>% unlist %>% as.numeric,
          MaximumYear=lapply(dat[[i]], FUN=function(x)x$SeaTemperature$MaximumYear) %>% unlist %>% as.numeric,
          Mean=lapply(dat[[i]], FUN=function(x)x$SeaTemperature$Mean) %>% unlist %>% as.numeric,
          Minimum=lapply(dat[[i]], FUN=function(x)x$SeaTemperature$Minimum) %>% unlist %>% as.numeric,
          MinimumYear=lapply(dat[[i]], FUN=function(x)x$SeaTemperature$MinimumYear) %>% unlist %>% as.numeric,
          MinimumAnomaly=lapply(dat[[i]], FUN=function(x)x$SeaTemperature$MinimumAnomaly) %>% unlist %>% as.numeric,
          MaximumAnomaly=lapply(dat[[i]], FUN=function(x)x$SeaTemperature$MaximumAnomaly) %>% unlist %>% as.numeric,
          MeanAnomaly=lapply(dat[[i]], FUN=function(x)x$SeaTemperature$MeanAnomaly) %>% unlist %>% as.numeric
          ) %>% as.data.frame
  )
}

# Hualien Buoy
hb <- MonthlyTemperatureFun(18) %>% subset((DataMonth=="4"|DataMonth=="8") & MinimumYear == "2012" & MaximumYear == "2012")
# Chenggong
cg <- MonthlyTemperatureFun(8) %>% subset((DataMonth=="4"|DataMonth=="8") & MinimumYear == "2012" & MaximumYear == "2012")
# Taitung Buoy
tb <- MonthlyTemperatureFun(34) %>% subset((DataMonth=="4"|DataMonth=="8") & MinimumYear == "2012" & MaximumYear == "2012")

tmp <- rbind(hb, cg, tb)

Aggregate by microhabitat

ta2 <- aggregate(b, by = list(ta$Site, ta$Month), FUN=sum)
names(ta2)[1:2] <- c("Site", "Month")
b2 <- ta2[, -1:-2]
rownames(b2) <- with(ta2, paste(Site, Month, sep="_"))

Estimate Hill Numbers

q0 <- iNEXT3D(t(b2), diversity = "TD", q = 0, nboot=1000)
q1 <- iNEXT3D(t(b2), diversity = "TD", q = 1, nboot=1000)
q2 <- iNEXT3D(t(b2), diversity = "TD", q = 2, nboot=1000)
#save(list=c("q0", "q1", "q2"), file="../Data/HillNumbersTD.rda")
load("../Data/HillNumbersTD.rda")

h <- rbind(q0$iNextEst$size_based,
           q1$iNextEst$size_based,
           q2$iNextEst$size_based)
h <- cbind(h, ta2[match(h$Assemblage, rownames(b2)), 1:2])
h$Month <- factor(h$Month, levels=c("4", "8"))

Fig. X. Size-based diversity accumulation curves based on spe richness (left), effective number of typical species (middle) and effective number of dominant species (right). The blue lines indicate the interpolated (rarefied) and red lines indicate extrapolated parts of the accumulation curves based on 1000 permutations. Dotted symbols indicate the observed diversity values.

# Minimum sample coverage for each order q
min(ldply(lapply(splitBy(~Order.q+Assemblage, subset(h, Order.q==0)), FUN=function(x)x[dim(x)[1],]))$SC)
## [1] 0.9957578
min(ldply(lapply(splitBy(~Order.q+Assemblage, subset(h, Order.q==1)), FUN=function(x)x[dim(x)[1],]))$SC)
## [1] 0.9957578
min(ldply(lapply(splitBy(~Order.q+Assemblage, subset(h, Order.q==2)), FUN=function(x)x[dim(x)[1],]))$SC)
## [1] 0.9957578
# SC = 0.99
h1 <- subset(h, Order.q ==0 & SC <= 0.9957578)

p1 <- ggplot()+
  geom_ribbon(data=h1,
              aes(x=m, ymin=qD.LCL, ymax=qD.UCL, group=Assemblage, fill=Site), alpha=0.2)+
  geom_line(data= subset(h1, Method=="Rarefaction"), 
            aes(x=m, y=qD, group=Assemblage, colour=Site), size=0.8)+
  geom_line(data= subset(h1, Method=="Extrapolation"), 
            aes(x=m, y=qD, group=Assemblage, colour=Site), linetype=2, size=0.8)+
  geom_point(data= subset(h1, Method=="Observed"), 
             aes(x=m, y=qD, group=Assemblage))+
  scale_fill_viridis_d()+
  scale_color_viridis_d()+
  facet_wrap(~Month, scale="free")+
  labs(x = "Numbers of Individuals", y = expression(Species~Richness~(""^0*italic(D))),
       colour="Site")+
  theme_bw()%+replace% large %+replace% theme(axis.text.x = element_blank(), axis.title.x = element_blank())
h2 <- subset(h, Order.q ==1 & SC <= 0.9957578)

p2 <- ggplot()+
  geom_ribbon(data=h2,
              aes(x=m, ymin=qD.LCL, ymax=qD.UCL, group=Assemblage, fill=Site), alpha=0.2)+
  geom_line(data= subset(h2, Method=="Rarefaction"), 
            aes(x=m, y=qD, group=Assemblage, colour=Site), size=0.8)+
  geom_line(data= subset(h2, Method=="Extrapolation"), 
            aes(x=m, y=qD, group=Assemblage, colour=Site), linetype=2, size=0.8)+
  geom_point(data= subset(h2, Method=="Observed"), 
             aes(x=m, y=qD, group=Assemblage))+
  scale_fill_viridis_d()+
  scale_color_viridis_d()+
  facet_wrap(~Month, scale="free")+
  labs(x = "Numbers of Individuals", y = expression(Species~Richness~(""^1*italic(D))),
       colour="Site")+
  theme_bw()%+replace% large %+replace% theme(legend.position = "none")
p1/p2

p3 <- ggplot()+
  geom_ribbon(data=h1,
              aes(x=SC, ymin=qD.LCL, ymax=qD.UCL, group=Assemblage, fill=Site), alpha=0.2)+
  geom_line(data= subset(h1, Method=="Rarefaction"), 
            aes(x=SC, y=qD, group=Assemblage, colour=Site), size=0.8)+
  geom_line(data= subset(h1, Method=="Extrapolation"), 
            aes(x=SC, y=qD, group=Assemblage, colour=Site), linetype=2, size=0.8)+
  geom_point(data= subset(h1, Method=="Observed"), 
             aes(x=SC, y=qD, group=Assemblage))+
  scale_fill_viridis_d()+
  scale_color_viridis_d()+
  facet_wrap(~Month, scale="free")+
  labs(x = "Sample Completeness", y = expression(Species~Richness~(""^0*italic(D))),
       colour="Site")+
  theme_bw()%+replace% large %+replace% theme(axis.text.x = element_blank(), axis.title.x = element_blank())
p4 <- ggplot()+
  geom_ribbon(data=h2,
              aes(x=SC, ymin=qD.LCL, ymax=qD.UCL, group=Assemblage, fill=Site), alpha=0.2)+
  geom_line(data= subset(h2, Method=="Rarefaction"), 
            aes(x=SC, y=qD, group=Assemblage, colour=Site), size=0.8)+
  geom_line(data= subset(h2, Method=="Extrapolation"), 
            aes(x=SC, y=qD, group=Assemblage, colour=Site), linetype=2, size=0.8)+
  geom_point(data= subset(h2, Method=="Observed"), 
             aes(x=SC, y=qD, group=Assemblage))+
  scale_fill_viridis_d()+
  scale_color_viridis_d()+
  facet_wrap(~Month, scale="free")+
  labs(x = "Sample Completeness", y = expression(Species~Richness~(""^1*italic(D))),
       colour="Site")+
  theme_bw()%+replace% large %+replace% theme(legend.position = "none")
p3/p4

fr <- strsplit(q0$DataInfo$Assemblage, split="_") %>% ldply
names(fr) <- c("Site", "Month")
fr$Month <- factor(fr$Month, levels=c("4", "8"))
h1_out <- lapply(splitBy(~Assemblage, h1), FUN=function(x)tail(x, n=1)) %>% ldply
h2_out <- lapply(splitBy(~Assemblage, h2), FUN=function(x)tail(x, n=1)) %>% ldply
t0 <- cbind(cbind(fr, S.obs = q0$DataInfo$S.obs), tmp[c(1, 3, 5, 2, 4, 6),]) 
t0$Site <- factor(t0$Site, levels=c("A", "B", "C"))
t1 <- cbind(h1_out, tmp[c(1, 3, 5, 2, 4, 6),])
t2 <- cbind(h2_out, tmp[c(1, 3, 5, 2, 4, 6),])
h2 <- readRDS("../Data/H2.rds")
h2$Site <- factor(h2$Site, levels=c("A", "B", "C"))
pv <- round(summary(lm(S.obs~Mean, data=t0))$coefficients[2, 4], 3)
r2 <- round(summary(lm(S.obs~Mean, data=t0))$adj.r.squared, 2)
f <- round(summary(lm(S.obs~Mean, data=t0))$fstatistic[1], 2)
label = paste0("italic(F)['1,4'] ==", f, "*','~P==", pv, "*','~italic(R)^2 ==", r2)

p7 <- ggplot(data=t0, aes(x=Mean, y=S.obs, shape=Month, colour=Site))+
  geom_point(size=2, stroke=1.5)+
  stat_smooth(formula=y~x, method = "glm", aes(group=1), alpha=0.2, colour="gray", linetype=2)+
  annotate(geom = "text", x = 24, y = 12.5, hjust = 0, vjust = 1, label = label, parse = TRUE) +
  scale_color_viridis_d()+
  scale_shape_manual(values=c(1,19))+
  labs(x = expression(Temperature~(degree~C)), y = "Numbers of Species", title="a", shape="Month")+
  theme_bw() %+replace% large %+replace% theme(axis.text.x = element_blank(), axis.title.x = element_blank(), legend.position = "none")
pv <- round(summary(lm(qD~Mean, data=t1))$coefficients[2, 4], 3)
r2 <- round(summary(lm(qD~Mean, data=t1))$adj.r.squared, 2)
f <- round(summary(lm(qD~Mean, data=t1))$fstatistic[1], 2)
label = paste0("italic(F)['1,4'] ==", f, "*','~P==", pv, "*','~italic(R)^2 ==", r2)

p8 <- ggplot(data=t1, aes(x=Mean, y=qD, ymin=qD.LCL, ymax=qD.UCL, shape=Month, colour=Site))+
  geom_point(size=2, stroke=1.5)+
  geom_errorbar()+
  stat_smooth(formula=y~x, method = "glm", aes(group=1), alpha=0.2, colour="gray", linetype=2)+
  annotate(geom = "text", x = 24.5, y = 22, hjust = 0, vjust = 1, label = label, parse = TRUE) +
  scale_color_viridis_d()+
  scale_shape_manual(values=c(1,19))+
  labs(x = expression(Temperature~(degree~C)), y = expression(Species~Richness~(""^0*italic(D))), title="b", shape="Month")+
  theme_bw() %+replace% large %+replace% theme(axis.text.x = element_blank(), axis.title.x = element_blank())
pv <- round(summary(lm(qD~Mean, data=t2))$coefficients[2, 4], 3)
r2 <- round(summary(lm(qD~Mean, data=t2))$adj.r.squared, 2)
f <- round(summary(lm(qD~Mean, data=t2))$fstatistic[1], 2)
label = paste0("italic(F)['1,4'] ==", f, "*','~P==", pv, "*','~italic(R)^2 ==", r2)

p9 <- ggplot(data=t2, aes(x=Mean, y=qD, ymin=qD.LCL, ymax=qD.UCL, shape=Month, colour=Site))+
  geom_point(size=2, stroke=1.5)+
  geom_errorbar()+
  stat_smooth(formula=y~x, method = "glm", aes(group=1), alpha=0.2, colour="gray", linetype=2)+
  annotate(geom = "text", x = 24, y = 6, hjust = 0, vjust = 1, label = label, parse = TRUE) +
  scale_color_viridis_d()+
  scale_shape_manual(values=c(1,19))+
  labs(x = expression(Temperature~(degree~C)), y = expression(Shannon~Diversity~(""^1*italic(D))), title="c", shape="Month")+
  theme_bw()%+replace% large %+replace% theme(legend.position = "none")
pv <- round(summary(lm(H2~Mean, data=h2))$coefficients[2, 4], 3)
r2 <- round(summary(lm(H2~Mean, data=h2))$adj.r.squared, 2)
f <- round(summary(lm(H2~Mean, data=h2))$fstatistic[1], 2)
label = paste0("italic(F)['1,4'] ==", f, "*','~P==", pv, "*','~italic(R)^2 ==", r2)

p10<- ggplot(data=h2, aes(x=Mean, y=H2, ymin=X1, ymax=X2, shape=Month, colour=Site))+
  geom_point(size=2, stroke=1.5)+
  geom_errorbar()+
  stat_smooth(formula=y~x, method = "glm", aes(group=1), alpha=0.2, colour="gray")+
  annotate(geom = "text", x = 24.5, y = 0.95, hjust = 0, vjust = 1, label = label, parse = TRUE) +
  scale_color_viridis_d()+
  scale_shape_manual(values=c(1,19))+
  labs(x = expression(Temperature~(degree~C)), y = "H2' index", title="d", shape="Month")+
  theme_bw()%+replace% large %+replace% theme(legend.position = "none")
(p7+p8)/(p9+p10)

p11 <- ggplot(data=t0, aes(x=Site, y=S.obs, fill=Month))+
  geom_bar(stat = "identity", position=position_dodge())+
  scale_fill_viridis_d()+
  labs(x = "Site", y = "Numbers of Species", title="a", shape="Month")+
  theme_bw() %+replace% large %+replace% theme(axis.text.x = element_blank(), axis.title.x = element_blank(), legend.position = "none")

p12 <- ggplot(data=t1, aes(x=Site, y=qD, ymin=qD.LCL, ymax=qD.UCL, fill=Month))+
  geom_bar(stat = "identity", position=position_dodge())+
  geom_errorbar(stat = "identity", position=position_dodge(), linetype=2)+
  scale_fill_viridis_d()+
  labs(x = expression(Temperature~(degree~C)), y = expression(Species~Richness~(""^0*italic(D))), title="b", shape="Month")+
  theme_bw() %+replace% large %+replace% theme(axis.text.x = element_blank(), axis.title.x = element_blank())

p13 <- ggplot(data=t2, aes(x=Site, y=qD, ymin=qD.LCL, ymax=qD.UCL, fill=Month))+
  geom_bar(stat = "identity", position=position_dodge())+
  geom_errorbar(stat = "identity", position=position_dodge(), linetype=2)+
  scale_fill_viridis_d()+
  labs(x = "Site", y = expression(Shannon~Diversity~(""^1*italic(D))), title="c", shape="Month")+
  theme_bw()%+replace% large %+replace% theme(legend.position = "none")

p14<- ggplot(data=h2, aes(x=Site, y=H2, ymin=X1, ymax=X2, fill=Month))+
  geom_bar(stat = "identity", position=position_dodge())+
  geom_errorbar(stat = "identity", position=position_dodge(), linetype=2)+
  scale_fill_viridis_d()+
  labs(x = "Site", y = "H2' index", title="d", shape="Month")+
  theme_bw()%+replace% large %+replace% theme(legend.position = "none")
(p11+p12)/(p13+p14)

library(knitr)
kable(t0)
Site Month S.obs StationNameEN StationLatitude StationLongitude DataMonth Maximum MaximumYear Mean Minimum MinimumYear MinimumAnomaly MaximumAnomaly MeanAnomaly
117 A 4 7 Hualien Buoy 24.03 121.63 4 28.0 2012 25.5 22.9 2012 -1.8 3.3 0.8
64 B 4 9 Chenggong 23.10 121.38 4 25.7 2012 24.0 22.7 2012 -2.6 0.4 -1.3
32 C 4 10 Taitung Buoy 22.72 121.14 4 29.0 2012 26.6 24.9 2012 -1.0 3.1 0.7
121 A 8 5 Hualien Buoy 24.03 121.63 8 31.1 2012 29.2 25.3 2012 -3.2 2.6 0.7
68 B 8 11 Chenggong 23.10 121.38 8 29.5 2012 28.1 26.9 2012 -1.6 1.0 -0.4
36 C 8 8 Taitung Buoy 22.72 121.14 8 30.0 2012 27.9 25.5 2012 -2.9 1.6 -0.5
kable(t1)
.id Assemblage m Method Order.q qD qD.LCL qD.UCL SC SC.LCL SC.UCL Site Month StationNameEN StationLatitude StationLongitude DataMonth Maximum MaximumYear Mean Minimum MinimumYear MinimumAnomaly MaximumAnomaly MeanAnomaly
117 A_4 A_4 211 Rarefaction 0 6.365286 5.627048 7.103523 0.9956038 0.9926322 0.9985754 A 4 Hualien Buoy 24.03 121.63 4 28.0 2012 25.5 22.9 2012 -1.8 3.3 0.8
64 B_4 B_4 298 Extrapolation 0 9.517606 7.651172 11.384039 0.9956365 0.9896097 1.0000000 B 4 Chenggong 23.10 121.38 4 25.7 2012 24.0 22.7 2012 -2.6 0.4 -1.3
32 C_4 C_4 444 Rarefaction 0 7.611364 6.480387 8.742340 0.9954990 0.9936961 0.9973019 C 4 Taitung Buoy 22.72 121.14 4 29.0 2012 26.6 24.9 2012 -1.0 3.1 0.7
121 A_8 A_8 138 Extrapolation 0 5.282364 3.289271 7.275457 0.9956636 0.9855805 1.0000000 A 8 Hualien Buoy 24.03 121.63 8 31.1 2012 29.2 25.3 2012 -3.2 2.6 0.7
68 B_8 B_8 1428 Extrapolation 0 14.932366 9.506878 20.357854 0.9957578 0.9924121 0.9991035 B 8 Chenggong 23.10 121.38 8 29.5 2012 28.1 26.9 2012 -1.6 1.0 -0.4
36 C_8 C_8 352 Rarefaction 0 6.880689 5.686621 8.074758 0.9953146 0.9929312 0.9976980 C 8 Taitung Buoy 22.72 121.14 8 30.0 2012 27.9 25.5 2012 -2.9 1.6 -0.5
kable(t2)
.id Assemblage m Method Order.q qD qD.LCL qD.UCL SC SC.LCL SC.UCL Site Month StationNameEN StationLatitude StationLongitude DataMonth Maximum MaximumYear Mean Minimum MinimumYear MinimumAnomaly MaximumAnomaly MeanAnomaly
117 A_4 A_4 211 Rarefaction 1 1.883177 1.681757 2.084597 0.9956038 0.9925190 0.9986886 A 4 Hualien Buoy 24.03 121.63 4 28.0 2012 25.5 22.9 2012 -1.8 3.3 0.8
64 B_4 B_4 298 Extrapolation 1 5.129220 4.604312 5.654129 0.9956365 0.9896115 1.0000000 B 4 Chenggong 23.10 121.38 4 25.7 2012 24.0 22.7 2012 -2.6 0.4 -1.3
32 C_4 C_4 444 Rarefaction 1 1.704697 1.604944 1.804450 0.9954990 0.9936266 0.9973714 C 4 Taitung Buoy 22.72 121.14 4 29.0 2012 26.6 24.9 2012 -1.0 3.1 0.7
121 A_8 A_8 138 Extrapolation 1 2.936496 2.495292 3.377701 0.9956636 0.9859488 1.0000000 A 8 Hualien Buoy 24.03 121.63 8 31.1 2012 29.2 25.3 2012 -3.2 2.6 0.7
68 B_8 B_8 1428 Extrapolation 1 3.509970 3.259010 3.760930 0.9957578 0.9921777 0.9993379 B 8 Chenggong 23.10 121.38 8 29.5 2012 28.1 26.9 2012 -1.6 1.0 -0.4
36 C_8 C_8 352 Rarefaction 1 3.880164 3.705037 4.055290 0.9953146 0.9928583 0.9977709 C 8 Taitung Buoy 22.72 121.14 8 30.0 2012 27.9 25.5 2012 -2.9 1.6 -0.5